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Abstract 

The homogeneous and heterogeneous nucleation of a Lennard-Jones Hquid is investigated using 
the umbrella sampling method. The free energy cost of forming a nucleating droplet is determined 
as a function of the quench depth, and the saddle point nature of the droplets is verified using 
an intervention technique. The structure and symmetry of the nucleating droplets is found for 
a range of temperatures. We find that for deep quenches the nucleating droplets become more 
anisotropic and diffuse with no well defined core or surface. The environment of the nucleating 
droplets form randomly stacked hexagonal planes. This behavior is consistent with a spinodal 
nucleation interpretation. We also find that the free energy barrier for heterogeneous nucleation is 
a minimum when the lattice spacing of the impurity equals the lattice spacing of the equilibrium 
crystalline phase. If the lattice spacing of the impurity is different, the crystal grows into the bulk 
instead of wetting the impurity. 



I. INTRODUCTION 



Although nucleation from a supercooled liquid has been the subject of extensive sim- 
ulations [D El El 11 El E], theory [3 El El ID], and experiments [III [12], the nature of 
the nucleating droplet in supercooled liquids is not well understood, especially for deep 
quenches. For shallow quenches (near coexistence), classical nucleation theory applies. For 
deeper quenches, nucleation is affected by the proximity to the liquid-solid spinodal for 
systems with long-range interactions [10]. The spinodal represents the limit of stability of 
the metastable liquid and is well defined only in the limit of an infinite interaction range. 
However, spinodal-like effects have been found for deep quenches in systems with inter- 
mediate and short-range interactions [H [131 [S] • Spinodal nucleation theory predicts that 
the decrease of the surface tension of the droplets as the spinodal is approached makes the 
nucleating droplets diffuse and fractal-like. Moreover, the symmetry of the nucleating drop- 
lets is not necessarily the same as the symmetry of the stable phase, and the symmetry of 
the nucleating droplets in three dimensions is either body-centered cubic (bcc) or randomly 
stacked hexagonal planes. Trudu et al. [H] studied nucleation of a Lennard- Jones liquid us- 
ing transitional path sampling, and found a crossover from classical to spinodal-like behavior 
for deeper quenches. In particular, they observed that the nucleating droplets become less 
compact and spherical, but did not analyze the structure and symmetry of the nucleating 
droplets. 

The study of heterogeneous nucleation, that is, nucleation that occurs on impurities, 
is of much practical importance because most nucleation events that occur in nature are 
heterogeneous. Examples include nucleation on a container wall [L5l and nucleation of 
proteins in porous media [16] . Many experiments [T71 [18] and simulations [HI [20l [211 [22] have 
been done to study heterogeneous nucleation. Existing theories of heterogeneous nucleation 
are mostly phenomenological and natural extensions of classical nucleation theory [231 [2H 
[25] . It is known that the presence of impurities can lower the free energy barrier of nucleation 
by as much as several orders of magnitude [H] . The effectiveness of an impurity to decrease 
the nucleation barrier is determined by properties such as the shape of the impurity and 
the surface tension between the substrate and the metastable liquid. Page and Sear [20] 
studied heterogeneous nucleation in porous media using the Ising model and found that 
a pore which is approximately the size of the critical nucleus is optimal for decreasing 
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the nucleation barrier. Heterogeneous nucleation on a structureless solid surface has also 
been simulated [T3]. However, the effects of the microscopic properties of the impurities on 
nucleation have not been well characterized. 

In this paper, we will study homogeneous and heterogeneous nucleation in supercooled 
Lennard- Jones liquids using the umbrella sampling method. For homogeneous nucleation, 
we find spinodal effects for deep quenches by analyzing the structure of the nucleating drop- 
lets. In particular, the nucleating droplets are found to become more anisotropic and diffuse 
with no well defined core or surface. The droplets and their immediate environment form 
randomly stacked hexagonal planes, which is consistent with the spinodal nucleation picture. 

To study heterogeneous nucleation a fixed impurity consisting of particles that form a 
hexagonal plane is added to the system. We find that the impurity whose lattice spacing is 
equal to the lattice spacing of the equilibrium crystalline phase is most effective in lowering 
the free energy barrier of nucleation. We also find that when the lattice spacing of the 
impurity is different than the optimal spacing, the crystal prefers to nucleate on the newly 
formed crystal (grow into the bulk) instead of wetting the impurity. 

This paper is organized as follows. Section [TT] describes the simulation details and the 
intervention technique which we use to test the saddle point nature of the nucleating droplets. 



Section [nT] presents the simulation results of homogeneous nucleation, and Sec. [IV| discusses 
our results on heterogeneous nucleation. 



II. SIMULATION DETAILS 



The three-dimensional system of interest consists of = 4000 particles with density 
p = 0.95 interacting via the Lennard- Jones potential. Periodic boundary conditions are 
used. We adopt dimensionless units so that lengths and energies are given in terms of the 
Lennard- Jones parameters a and e. We first prepared a liquid at T = 1.20, which is above the 
coexistence temperature (Tm ~ 1.15 [26]), by melting a perfect fee crystal; this simulation 
is done using the Metropolis algorithm at constant volume. The system is equilibrated for 
50,000 Monte Carlo steps per particle (mcs) in the liquid phase before the quench. The 
system is then quenched by rescaling the temperature by a factor of 0.999 every 20 mcs. 

Because the probability of nucleation is very small, we used the umbrella sampling 
method [271 EHl EHl 1301 EI]. We denote (p as the order parameter, which we choose to 



3 



be the number of particles in the largest (solid-like) cluster. (The definition of the solid-like 
clusters is discussed in Sec. II A[ ) The free energy G(0) is calculated from the relation 



G(0) = -fcBTlnP(0), (1) 

where -P(0) is the probability density of 0. 

In the umbrella sampling method, the system is sampled according to the total energy 
V = + Vb(0), where V is the original potential energy of the system and Vb(0) is the 
bias potential. The probability distribution P(0) is sampled according to the total density 
operator p = e~^^ e~^^'°^'^^ = poe~^^^^'^\ which is the product of the original density operator 
Po and the weight function due to the bias e~^^^^'^\ The original distribution P(0) can be 
determined by 

P(0) = P(0)e^^''('^). (2) 
Hence, the free energy G{4>) in Eq. (|g can be calculated by 

G(0) = -A;BTlnP(0)-H(0). (3) 
As in Ref. [271 the potential bias has the form 

VM = lki<P-hy- (4) 

The constant k = 0.05 determines the width of the sampling window and yields A(f) ^ 15. 
We consider a sequence of values of 0o starting from size and increasing by ten particles, 
that is, (pQ = 0, 10, 20, .... Because the width of the window A(f) ^ 15, the choice of ten 
particles means that the sampling windows overlap. Before collecting data for the probability 
P(0) for each value of (po, the system is equilibrated for 10,000 mcs. The values of are then 
sampled for 100,000 mcs. To save equilibration time a configuration for the current value of 
00 is used as the initial condition for the next value of 0o. Because determining the size of 
the largest cluster is computationally expensive, we make trial moves of 5 mcs using only 
the Lennard- Jones potential without the bias potential (|4]), and then accept or reject these 
trial moves using only the bias potential in Eq. (|4]). 

A. Cluster analysis 

Unlike Ising/Potts models [32] there is no rigorous definition of clusters in a continuous 
particle system. Instead we are forced to rely on our intuition to identify the solid-like 



particles. We use the local bond-order analysis introduced by Steinhardt et al. [33] and 
developed by Frenkel and co-workers ^T\. We define the {21 + 1) component complex vector 
Qimii) for particle i: 

ni 

Qlmi'^ = —^yimirij), (5) 
' 0=1 

where the sum is over the rii nearest neighbors of particle i and l^m(^ij) is the spherical 
harmonic as a function of the unit direction vector r^j between particle i and its rii neighbors. 
The nearest neighbors of a given particle are defined to be within the distance 1.4, which 
corresponds to the position of the first minimum of the radial distribution function g{r) of 
the crystalline phase at the same density and temperature. 

It has been shown that Z = 6 is a good choice for characterizing the structures of crys- 
tals [28j. The rotational invariants Wi{i), Wq{i), qi{i), and q^ii) are defined as 

E \-qi^m"\ (6) 



An 
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m=—l 



and 



with 



Wl{l) 



= 7^1 ,_ (7) 
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mi+m2+m3=0 \ "^1 "^2 '^S / 

The quantity in brackets in Eq. ([s]) is the Wigner's 3j symbol. 

To define a solid-like particle, we first introduce the normalized quantity qim{i) as 

~ / -A _ Qlm{^) /q\ 

and form the dot product 

6 

^ij = Q6mii)q6mij)\ (10) 

m=—6 

where * indicates the complex conjugate. The dot product Cu = 1 by construction. Particles 
i and j are said to be coherent if the real part of the dot product Cij is greater than 0.5. A 
particle is considered to be solid-like if its number of coherent neighbors is greater than or 



equal to nb. We will chose nb = 11 for reasons that we will discuss in Sec. |11 B[ After finding 
all the solid-like particles, we identify the clusters using the criterion that any two solid-like 
particles that are nearest neighbors belong to the same cluster. 
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fraction 


0.75 


0.65 


0.60 


0.40 


0.70 


0.60 


0.55 


0.45 


0.65 


0.35 


0.53 


0.55 



TABLE I: Fraction of the copies of the nucleating droplets that grow after intervention (out of 
20). These results indicate that the clusters that are generated at the maximum of the free energy 
correspond to nucleating droplets. 

As discussed in Ref. [31], the distribution of the invariants g4(i), WQ{i), and qQ{i) can be 
used to characterize the symmetry of a group of particles, for example, a cluster. For a given 
group of particles we first determine the histogram of each invariant. The three histograms 
are then rescaled so that they do not overlap and are combined to form the histogram 
H. The histogram H corresponding to a group of particles is decomposed in terms of the 
histograms of the invariants corresponding to each symmetry of interest, namely, 

H = ficcHicc + /bcc-f^bcc + /hcp-f^hcp- (H) 

For example, H^^c is determined for a system with fee symmetry in which a certain amount 
of randomness has been introduced [3l]. The coefficients corresponding to each symmetry, 
/fee, /bee, and /hep, are found by minimizing the quantity 

{H — /fcc-f^fee — /bee-f^bec — /hep-f^hep)^, (12) 

with the constraints that /fcc, /bee, /hep > and /^c + /bee + /hep = 1- The coefficients 
/fee, /bee and /hep are indications of the composition of each symmetry associated with the 
group particles. 

B. The intervention method 

Because the cluster analysis involves parameters such as the minimum number of coherent 
neighbors, the computed size of the clusters depends somewhat on the values chosen for the 
parameters. To determine if our choices are consistent, we use the fact that the nuclea- 
ting droplets correspond to the maximum of the free energy barrier and should be saddle 
point objects; that is, the nucleating droplets should grow or shrink with approximately 50% 
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FIG. 1: (Color online) The size of the largest cluster as a function of the time after intervention 
for four trials. After iobs = 10000 mcs, the largest droplet in each trial is determined to grow or 
shrink based on its size and the position of its center of mass compared to the largest cluster in 
the original configuration. The size of the nucleating droplet in this case is ~ 80 particles (for 



probability under a small perturbation. In practice, the parameters in the cluster analysis 
are first chosen and the umbrella sampling procedure is performed. Then the intervention 
method is used to test if intervention causes the clusters corresponding to the maximum 
of the free energy to grow or shrink with ^ 50% probability. If these clusters do not, the 
parameters in the cluster analysis are modified until approximately 50% growth probability 
is achieved for the modified clusters at the new free energy maximum. 

The intervention method we adopt is similar to what has been used to study nucleation in 
the Ising model [35]. To test if a cluster is a nucleating droplet, we stop the simulation and 
make many copies of the system. Each copy is restarted with a different random number 
seed without the potential bias. We then determine if the largest cluster in each copy grows 
at approximately the same place at approximately the same time as the original. After the 
time tobs, both the size and location of the cluster are examined. If the size of the cluster is 
larger than its original size and the center of mass is within a distance r* from the cluster 
in the original configuration, the cluster is said to grow. The role of r* is to ensure that the 
cluster is the same as the one in the original configuration. We choose tobs = 10000 mcs. The 
distance r* should correspond to the size of the cluster; typically 50% of the linear spatial 



T = 0.65). 
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extent of the cluster is sufficient to decide whether the cluster kept its identity during tobs- 
Because intervention is very time consuming, we made 20 copies of each configuration. The 
evolution of the size of the largest cluster after intervention for four trials is plotted in Fig. [T] 
for T = 0.65. Note that the largest cluster has grown in some trials and shrunk in others. By 
determining the frequency of the successful trials we can estimate the probability of growth 
of a cluster. The fraction of trials (out of 20 trials) for which the cluster grows is listed in 
Table [T] for different temperatures. We found that setting ric = 11 gives consistent results. 



III. HOMOGENEOUS NUCLEATION 



We determined G{(j)) at several temperatures between T = 0.75 and T = 0.53. As 



expected, G{(j)) exhibits a minimum and a maximum (see Fig. 2(a)). The minimum of G 
corresponds to the metastable supercooled liquid, where solid-like particles are present due to 
thermal fluctuations. The maximum occurs at the size of the nucleating droplet (recall that 
we have chosen the order parameter to be the size of the largest cluster). The free energy 
difference of the maximum and the minimum can be interpreted the free energy barrier 
AG for nucleation [27J. We find that the nucleation barrier decreases with temperature 



(see Fig. 2(b) ) and vanishes at T ^ 0.53. This vanishing of the free energy barrier raises 
questions on whether or not it corresponds to a spinodal. 

The spinodal is usually defined as the sharp boundary between the metastable and un- 
stable states [37]. In particular, the spinodal is a thermodynamic transition that acts as 
a line of critical points. (In Ising models the spinodal corresponds to a divergent isother- 
mal susceptibility [38].) We will refer to this interpretation as the classical spinodal. The 
spinodal as defined in this way is present only in mean-field systems such as those with 
infinite range interactions. Systems with large but finite range interactions can exhibit 
pseudospinodal effects [39]. To test if T = 0.53 corresponds to the classical spinodal we also 
simulated the system at this temperature by a standard Metropolis algorithm. By tracking 
the size of the largest cluster we found that the lifetime of the metastable state is in the 
range [1 x 10^, 6 x 10^] mcs, which implies that the free energy barrier to nucleation has 
not vanished at T = 0.53. Hence the vanishing of the free energy barrier found by umbrella 
sampling does not necessarily correspond to a classical spinodal. Moreover, we found (see 
Table |l]) that the droplets found by umbrella sampling appear to be saddle point objects as 
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(a) The Gibbs free energy. 



(b) The nucleation barrier. 



FIG. 2: (a) The Gibbs free energy as a function of the order parameter (p^ the size of the largest 
cluster in the system, at T = 0.55. The minimum of G corresponds to the metastable supercooled 
liquid phase, and the maximum corresponds to the nucleating droplet. The free energy minimum 
is at (/) « 7, showing that there are some solid-like particles in the supercooled liquid state, (b) 
The nucleation barrier as a function of the temperature. The barrier vanishes at T ~ 0.53. 

determined by the intervention method (without the bias potential). A possible explanation 
is that the interpretation of the umbrella sampling results for P(0) assumes that clusters 
whose size are comparable to the nucleating droplet are rare [21] • This assumption is not 
applicable for T ^ 0.53 because there is typically more than one large cluster of comparable 
size in the system. We will investigate this assumption and other possible explanations in 
future work. 

Figure [3] shows snapshots of the nucleating droplet at different temperatures. Note that 
the droplets are compact for moderate supercooling and become more diffuse for deeper 
quenches. This qualitative observation is consistent with Ref. [131 We will analyze the 
structure of the nuclei in the following. 

A. The structure of the nucleating droplets 

To measure the compactness of a nucleating droplet, we determine its density profile 
p(r), which is defined in terms of the mean number of particles N{r) in the spherical shell 
between r and r + dr 



N{r) = p{r)4:7rr'^dr. 



(13) 
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(c) T = 0.60. (d) T = 0.55. 

FIG. 3: (Color online) Snapshots of the nucleating droplet at different temperatures. 



Here r is measured from the center of mass of a nucleating droplet. 

Figure |4] shows p(r) averaged over 1000 nucleating droplets at the value of 0o correspond- 
ing to the maximum of G{(j)). For T = 0.75, p(r) has a plateau for small r, meaning that the 
droplet has a well defined core. The decrease of p(r) for larger r indicates that there is an 
interface between the core and the liquid environment. At T = 0.55, the plateau disappears, 
and the density slowly decreases from the core to the surface, indicating that the nucleating 
droplet is diffuse. 

To quantify the anisotropy of the nucleating droplets, we calculate the moment of inertia 
tensor associated with each droplet, 

n 

lap = ^{riSaf3 - ri,ari,(i), (14) 
i=l 

where rf = "^^ri^ari^a, i labels the particles, and a and (3 label the components of r. The 
square root of the eigenvalues of I^p define the principal radii of the ellipsoid characterizing 
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FIG. 4: The density profile p{r) of tlie nucleating droplets for various temperatures. At T = 0.75, 
the nucleating droplet has a well defined core corresponding to the plateau of p{r). Note that the 
density at the core is slightly higher than 0.95, the mean density. The decrease of p(r) for larger 
r implies that there is a well defined interface between the core and the liquid environment. At 
lower temperatures, the plateau disappears and the density changes gradually from the core to the 
surface, suggesting that the nucleating droplet are more diffuse. 

the droplet. The orientation of each individual nucleating droplet is found to be random; 
that is, the long axis of the computed ellipsoid points in random directions independent of 
the orientations of the simulation cell. We can characterize each droplet 's anisotropy by 
calculating the ratio of the maximum and minimum principal radii (denoted by Amax and 
Ajnin respectively). This ratio is one for a perfectly spherical droplet and is greater than one 
if the droplet is anisotropic. Figure |5] shows the ratio (averaged over 1000 nucleating droplets 
at each temperature) as a function of the temperature. The ratio is close to one for shallow 
quenches, meaning that the nucleating droplets are close to spherical. The increase of the 
ratio at lower temperatures indicates that the nucleating droplets become more anisotropic. 
The anisotropic character of the nucleating droplets is important in the calculation of the 
nucleation barrier even in classical nucleation theory [Ti] . 

Figure [6] shows that the number of particles Uc in the nucleating droplet and the radius of 
gyration Rg decrease as the temperature is decreased. Both quantities are predicted to first 
decrease as the temperature is lowered from coexistence and then begin to increase as the 
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FIG. 5: The ratio of the maximum and minimum of the eigenvalues of the moment of inertia 
tensor of the nucleating droplets. The increase of the ratio at low temperatures indicates that the 
nucleating droplets become more anisotropic. 

spinodal is approached [ID]. In particular, if the spinodal interpretation is applicable, simple 
scaling arguments [H] suggest that ric (in three dimensions) and Rg (in all dimensions) should 
increase as the spinodal is approached if the system is sufficiently close to the spinodal so 
that the core of the nucleating droplet has disappeared. The fact that Uc and Rg do not 
increase rapidly at lower temperatures in our simulations might be due to the nonexistence 




(a) Mean number of particles. (b) Radius of gyration and semimajor axis. 

FIG. 6: (a) The mean number of particles in the nucleation as a function of T. (b) The radius of 
gyration Rg and mean semimajor axis Amax of the nucleating droplets. The data is averaged over 
1000 independent configurations. Note that the Amax increases as T is decreased below T « 0.65. 
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of spinodal effects and/or the underestimate of and Rg near tlie spinodal due to our ad 
hoc definition of sohd-hke particles. A more hkely explanation is that because the nucleating 
droplets near the spinodal are anisotropic and effectively two-dimensional, the simple scaling 



arguments do not apply. In Fig. 6(b) we plot the semimajor axis Amax of the nucleating drop- 
lets. Note that Amax does show the expected behavior. In addition, the scaling arguments [H] 
in two dimensions suggest that is either a constant or is logarithmically divergent as the 
spinodal is approached. More work is needed to understand the temperature-dependence of 
Rg, and Amax in the intermediate region where neither the classical nor spinodal picture 
is applicable. 



B. Symmetry of the nucleating droplets 

The symmetry of the nucleating droplets is analyzed using the method discussed in 



Sec. II A At each temperature, we obtain the parameters /fcc, /bcc, and /hep for the nuc- 
leating droplets averaged over 1000 independent configurations. The fitting parameters are 
listed in Table [TT] and are plotted in Fig. [7j As the temperature is decreased, the fcc com- 
ponent decreases and the bcc and hep components increase. The mixture of fcc and hep 
signifies the occurrence of the rhcp structure [12] and is consistent with the picture of stacked 
hexagonal planes [TU] . 

We also calculated the symmetry of the particles in a spherical shell between r and r + Ar, 
where r is measured from the center of mass of the cluster and Ar is the thickness of the shell 
(Ar = 0.2). Figure Is] shows the component of each structure as a function of r for various 
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/fee 


/hep 


/bcc 


0.75 


0.62 


0.38 


0.00 


0.70 


0.50 


0.48 


0.02 


0.65 


0.46 


0.52 


0.02 


0.60 


0.28 


0.66 


0.06 


0.55 


0.26 


0.68 


0.06 



TABLE II: Values of the coefficients corresponding to the symmetry of the nucleating droplets for 
different temperatures. 
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FIG. 7: The coefficient / of each symmetry of the nucleating droplet as a function of temperature. 
The results are averaged over 1000 independent configurations. As the temperature is decreased, 
the hep and bcc components increase and the fee component decreases. 

temperatures. At T = 0.75, the core of the nucleating droplet is mostly fee. Away from 
the center, the fee component decreases and the bcc and liquid component increases. At 
the surface, the bcc component levels off to /bcc ~ 0.1. The fact that the nucleating droplet 
is composed of an fee core and a bcc halo agrees with previous results [31j. At T = 0.55 
the nucleating droplet is mainly a mixture of fee and hep, with a slight increase of the bcc 
component. The increase of bcc symmetry and decreased distinction between the bulk and 
the surface for deep quenches is in agreement with the spinodal nucleation picture |10j . 

We also examined the structure of the particles in the nucleating droplet and its local 
environment, which consists of particles that are nearest neighbors of any particle in the 
droplet. Although the particles in the nucleating droplets by themselves do not seem to 
form a visually identifiable structure, the nucleating droplets and the surrounding particles 
in the liquid phase together form hexagonally stacked planes (see Fig. [o]). 




IV. HETEROGENEOUS NUCLEATION 



To study heterogeneous nucleation, an impurity of m x n Lennard- Jones particles in a 
hexagonal plane is placed into the system, where m and n are the number of particles in 



the X and y directions (see Fig. 10). The z direction is perpendicular to the plane of the 
impurity. The positions of the particles in the impurity are fixed during the simulation. The 
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(c) T = 0.60. 



(d) T = 0.55. 



FIG. 8: The component of each symmetry as a function of the distance r from the center of 
mass of the nucleating droplets at different temperatures. The results are averaged over 1000 
configurations. At T = 0.75, the core of the nucleus is predominately fee, and the surface shows 
some bcc and liquid symmetry. At T = 0.55, the fee symmetry is less dominant for small r, and 
the bcc component for small r is close to its value for larger r. 



impurity is characterized by its lattice spacing a and total area A = ^mna^. The efficiency 
of the impurity is measured by the height of the nucleation barrier AG, which we compute 
as before using the umbrella sampling method. 

If the system is crystallized homogeneously after a quench to T < Tm, the position of the 
ffist peak of the radial distribution function g{r) is ~ 1.09. We take the value Og = 1.09 
as the lattice spacing of the solid phase. The temperature is quenched to T = 0.75, which 
corresponds to the region where classical nucleation applies in the absence of an impurity (see 



Sec. III). At T = 0.75 the free energy barrier of homogeneous nucleation is AG/k-^T ^ 40 
in the absence of an impurity with Uc ~ 300. In all of our simulations of heterogeneous 
nucleation, nucleation always occurs on the impurity if it is present. 
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(a) (b) 

FIG. 9: (Color online) (a) Snapshot of a nucleating droplet (dark) and particles in the local 
environment (light) at T = 0.55. (b) The layer of particles corresponding to the plane indicated 
by the dashed line in (a) form a hexagonal structure. 



FIG. 10: Sketch of the impurity used to study heterogeneous nucleation. The area of the impurity 
is yl = ^mna^, where m and n are the number of particles in the x and y directions respectively. 



We studied the dependence of the nucleation barrier on the lattice spacing for fixed area 
A, which is chosen to be A = ^(25 x 1.09^), the area of a 5 x 5 impurity with lattice spacing 
Os = 1.09. The impurities have different values of m and n so that their lattice spacings are 
less than and greater than as (see Table III). The free energy barrier is lowest at a ~ (see 



Fig. 11), that is, an impurity is most efficient in lowering the nucleation barrier if its lattice 



spacing is the same as that of the crystalline phase. 

By investigating snapshots of the nucleating droplets we found that in general the nuc- 
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m X n 


lattice spacing a 


AG/kBT 


6x6 


0.908 


37 


5x6 


0.995 


22 


5x5 


1.090 


3.3 


4x5 


1.218 


20 


4x4 


1.360 


30 



TABLE III: Values of m and n used to study the effect of changing the lattice spacing. The area 
of the impurity is fixed. 
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FIG. 11: The free energy barrier as a function of the lattice spacing a for fixed area X 
1.09^). The nucleation barrier is a minimum at a = a^. 



leating droplet forms on both sides of the impurity for a ^ (see Fig. |12(a) ). For a 7^ Og, 
once the impurity initiates nucleation, it becomes preferential to add particles on the newly 
formed droplet than on the impurity, making the droplet grow into the bulk instead of 
wetting the other side of the impurity. In both cases, the nucleating droplets grow on 
the impurity by forming layered planes, with each layer being parallel to the plane of the 



impurity. Inside each layer, the solid-like particles form a hexagonal structure (Fig. 13). 
The lattice spacing of the nucleating droplets is approximately irrespective of the lattice 
spacing of the impurity. That is, once a nucleus is formed, its lattice spacing becomes very 
close to the optimal spacing and the droplets form on the newly formed layer rather than 
on the impurity. 

A simple way to characterize how the nucleating droplet wets the impurity is to measure 
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(a) a = 1.09. (b) a = 0.908. 

FIG. 12: (Color online) Snapshots of nucleating droplets in the presence of an impurity. The dark 
particles are the impurity; the size of the nucleating droplets is 146 and 250 particles respectively. 
For a = 1.09, the nucleating droplet grows on both sides of the impurity. For a = 0.908, once the 
impurity initiates nucleation it is preferential to add particles on the newly formed droplet than 
wetting the other side of the impurity. 




(a) Impurity and first layer of droplet. (b) First layer of droplet only. 

FIG. 13: (Color online) Snapshots of the first layer of the nucleating droplet in the presence of 
a 5 X 5 impurity at a = Og. The layer of the nuclei lies on top of the impurity (dark particles). 
Within the layer the solid-like particles form a hexagonal structure, (b) Same snapshot as (a) with 
particles in the impurity not shown. 
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FIG. 14: Number of particles in the nuclei as a function of z. The solid line is for a 5 x 5 impurity 
with a = 1.09 and the dotted line is for a 6 x 6 impurity with a = 0.908. 



its profile. In Fig. 14 we show the number of particles in the nucleating droplet n{z) as a 
function of z, where z is the distance of a particle from the plane of the impurity. The sharp 
peaks of n{z) indicates that the nuclei form a layered structure. For a = 1.09 n{z) gradually 
drops to zero as z increases. In contrast, n{z) is almost fiat for a = 0.908. The protrusion of 
the nuclei into the bulk for a = 0.908 is consistent with the fact that the nucleating droplet 
wets the impurity at a = and grows into the bulk when a ^ ag. 

Contour plots of the profile were obtained by projecting the density of the nuclei onto 



the x-z plane. Figure [15] shows contour plots for impurities with a = 1.09 and a = 0.908 
respectively. The smaller contact angle for the impurity with a = 1.09 is consistent with the 
fact that it is preferential to nucleate on the impurity when its lattice spacing is optimal. 



V. CONCLUSIONS 



We have studied the homogeneous and heterogeneous nucleation of Lennard- Jones liquids 
using the umbrella sampling method. By analyzing the symmetries of the nucleating drop- 
lets, we found that for deep quenches the nucleating droplets are more diffuse and anisotropic 
with no well defined core or surface; the nucleating droplets and the corresponding liquid 
environment form randomly stacked hexagonal planes. These results are consistent with the 
spinodal nucleation picture. For heterogeneous nucleation, we found that the droplets grow 
on an impurity of hexagonal plane by layers, and the solid-like particles in each layer form 
a hexagonal structure. For fixed area of the impurity, the free energy barrier of nucleation 
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FIG. 15: (Color online) Contour plots of the nucleating droplets in the presence of an impurity. 
The impurity is at z = — 1 and the density of the nuclei is projected onto the x-z plane. The plots 
are obtained by averaging over 200 independent configurations of the nucleating droplets, (a) The 
lattice spacing of the impurity is a = 1.09, with ric = 150. (b) The The lattice spacing is a = 0.908 
with Uc = 300. The smaller contact angle of the impurity with a = 1.09 is consistent with the fact 
that it is preferential to nucleate on the impurity when its lattice spacing is optimal. 

is a minimum when the lattice spacing of the impurity is equal to Og, the lattice spacing of 
the equilibrium crystalline phase. The lattice spacing of the nuclei is equal to even when 
the lattice spacing of the impurity is different than a^, and it is favorable for the nucleating 
droplets to grow into the bulk instead of wetting the impurity. 
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